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ABSTRACT 


For the spacecraft and aircraft designers, the ability to change and control the shape of 
the structures has been a challenging problem. For a spacecraft, it is highly desirable to 
change the reflector shape in orbit to conqwnsate for the sur&ce errors and also to 
perform antenna beam shaping . In the case of aircraft, the sh£^ control of propellers, 
helicopter rotor blades, and aircraft wings can increase efficiency and maneuverability 
and reduce vibration and noise of the vehicles. In the present work, the sht^ control of 
fiber reinforced composite plates with embedded piezoelectric actuators is investigated. 

In the present study, a finite element formulation is developed for modeling a 
laminated composite plate that has distributed piezoelectric actuators and sensors 
subjected to both mechanical and electrical loads. A single, higher order, shear 
deformation theory with Hamilton's principle is used to formulate the equations of 
motion. The model represents the parabolic distribution of transverse shear stresses and 
the non-linearity of in-plane displacements across the thickness. The model is valid for 
both segmented and continuous piezoelectric elements, which can be either surface 
bonded or embedded in the laminated plate. A four-node, bilinear, isoparametric, 
rectangular element with seven degrees of fi'eedom at each node is developed. The 
electric potential is treated as a generalized electric coordinate like the generalized 
di^lacement coordinates at the mid-plane of the actuator layers. 

For shape control, an optimization algorithm, based on a finite element techniques, is 
presented for an optimal jq)plied voltage to each actuator to minimize the error between 
the desired shape and the actual shape. The error (objective) function is the mean square 
of the error between the point in the actual surftice and the corresponding point in the 
desired surfece. Based on these techniques, two computer programs have been developed, 
a finite element modeling of a composite plate with piezoelectric actuators and an 
optimization model of the actuator voltages for shape control. The present work 
demonstrates the feasibility of the application of the piezoelectric actuators for the shape 
control of con^site plates used in aero^ace structures. 
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I. INTRODUCTION 


Advanced structures with integrated self-monitoring and control capabilities are 
becoming incre asing ly iitqx)rtant due to the rapid development of ‘intelligent’ space 
structures. Those smart structures that have actuators distributed throughout are defined 
as adaptive, or actuated, structures. Exan 5 )les of such adaptive structures are conventional 
aircraft wings with articulated leading-and trailing-edge control sur&ces, robotic systems 
with articulated manipulators and end effectors and spacecraft antenna reflectors. 
Structures that have sensors distributed throughout are a subset referred to as sensory. 
These structures have sensors which can detect displacements, strains or other mechanical 
states or properties like electromagnetic states or properties, ten^)erature, heat flow, or 
damage. Applications of this technology might include damage detection in long life 
structures, or embedded or conformal RF antennas within a structure. The overlap 
structures, winch contain both actuators and sensors implicitly linked closed-loop 
control, are referred to as controlled structures. A subset of controlled structures are 
active structures, distinguished fi“om controlled structures by extensively distributed 
actuators, which have structural functionality and are part of the load bearing system. 

Intelligent structures are a subset of active structures that have highly distributed 
actuator and sensor systems with structural functionality and, in addition, distributed 
control functions and computing architecture [Ref 1]. 

Smart structures have applications due to their ability to change shape. Some of these 
are for spacecraft antennas, to compensate for surfiice error and thermal distortion to 
inprove antenna performance, and to change the antenna beam shape in orbit. They may 
be used for and submarine and helicopter shape control as weU. An aeroelastic application 
to aircraft structures is quasi-static control of camber, dynamic control, and flutter 
suppression. Smart structures can be used in acoustic control by developing an adaptive 
structure in which the structural response can be modified with var 5 Tng input disturbances. 
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The piezoelectric mat erial is one of the smart materials that can be used as an 
actuator and sensor. Piezoelectricity refers to the phenomenon of generating an electric 
charge in a material when subjecting it to mechanical stress, which is known as a direct 
piezoelectric effect. The converse piezoelectric effect is described as an induced strain in 
response to an applied electric field. 

Piezoelectric properties occur naturally in some crystalline materials and can be 
induced in other polycrystalline materials. The distortion of the crystal domain produces 
the piezoelectric effect. The domain may be aligned/poled by the q>plication of a large 
electric field, usualfy at high ten: 5 )erature. Subsequent ^plication of the electric field will 
produce additive strains in the local domain, >\iiich translate into a global strain in the 
material. 

The piezoelectric effect was discovered in 1880 by Pierre and Jacques Curie. The 
direct piezoelectric effect has been used for a long time in sensors such as accelerometers. 
Use of the converse effect, however, has imtfl recently been restricted to ultrasonic 
transducers. Barium titanate, discovered in 1940, was the first widely used piezoceramic. 
Lead zirconate titanate (PZT) [Ref 2], discovered ini954, has now largely superseded 
barium titanate because of its stronger piezoelectric effect. 

A P-phase polymeric piezoelectric, polyvinylidene fluoride (PVDF), was initially 
discovered by Kawai in 1969 [Ref 3]. Raw polymeric PVDF (a-phase) is an electrical 
insulator, and it does not have any intrinsic piezoelectric properties. If the raw material 
is polarized during the manufecturing process, PVDF transforms to a P-phase , a tough 
and flexible semi-crystalline material, and can be made to strain either in one or two 
directions in the film plane. Since P-phase PVDF possesses a strong direct piezoelectric 
effect, it has been used in many transducer applications; e.g., sonar, medical ultrasonic 
equipment, robot tractile sensors, acoustic pick-ups, force and strains gages, etc. Due to 
its distinct characteristics, such as flexibility, disability, manufecturability, etc., PVDF is 
commonly used in such structural systems. 
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As a first step in the following work, a finite element analysis of a graphite/epoxy, 
fiber reinforced, composite, laminated plate is developed, using a sin 5 )le higher-order 
shear deformation theory [Ref 4]. Multi-layered con^wshes have found wide use in many 
weight-sensitive structures, such as aircraft, spacecraft antennas, and missile structural 
conqwnents, where high strength-to-weight and stffiiess-to-weight ratios are required. A 
laminate is a multi-layered conyx>site made up of several individual layers ( lamin ae) in 
vvluch the fibers are oriented in a predetermined direction to provide eflBciently the 
required strength and stifibess parameters in each laminae. 

For analy sis and design of structural laminates, a classical plate theory (CPT) [Refe. 
5 ,6,7, & 8] has been used in t^ch it is assumed that normals to the mid-plane before 
deformation remain straight and normal to the mid-plane after deformation (classical 
Kirchhoff hypotheses). This theory under-predicts deflections and over-predicts natural 
fi'equencies and bucklii^ loads. These results are due to the neglect of transverse shear 
strains in the classical theory. These errors are even higher for a plate made of advanced 
conqKtshes like graphite epoxy and boron-epoxy, whose elastic modulus to shear 
modulus ratio are very large (Le., of the order of 25 to 40 instead of 2.6 for a typical 
isotropic materials). Tl^se high ratios render classical theories inadequate for the analysis 
of conqwshe plates. An adequate theory must account for transverse shear strains. 

Many plate theories have been proposed to incorporate the influence of the transverse 
shear strains. In one of these, the Reissner-Mindlin plate theory [Ref 9], transverse shear 
and rotary inertia effects are included, and it contains a displacement field which accounts 
for linear, or higher-order, variations of mid-plane displacements through the thickness, 
but on the other side the deviation increases with increasing mode numbers. A theory for 
laminated isotropic plates by Stavsky [Ref 10]. has been generalized to laminated 
anisotropic plates by Yang, Norris et. al. [Ref 11]. Whitney [Ref 12] has presented an 
approximate method to incorporate the influence of shear deformation on plate deflection, 
in a flexural vibration and buckling analysis. Elasticity solutions by Pagano and his 
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associates [Refe. 13, 14, 15, & 16] indicated the inadequacy of the classical laminate 
theory. These shear deformation theories do not satisfy the conditions of zero transverse 
shear stresses on the top and bottom surfeces of the plate, and they require a shear 
correction to the transverse shear stifGoess. Three-dimensional theories of laminates, in 
which each layer is treated as a homogeneous anisotropic medium [Refe. 17, 18, & 19] 
are intractable as the number of layers becomes moderately large. 

Different highe r-order laminated plate theories have been proposed, t\bich account 
for the transverse shear strains. Examples of such theories are, Whitney and Pagano [Ref 
20], Whitney et. al [Refe. 21,22], Lo et al. [Ref 23] and Nelson et. al. [Ref 24] In these 
hi ghe r-order theories, an additional dependent unknown is introduced into the theory with 
each additional power of the thickness coordinate. 

A single two-dimensional theory of plates that accurately describes the global 
behavior of laminated plates seems to be a compromise between accuracy and ease of 
analysis. A single, higher-order theory described by Reddy [Ref 25], is such a theory, 
as it is accounts not only for transverse shear strain but also for a parabolic variation of 
the transverse shear strains through the thickness. Consequently, there is no need to use 
shear correction coefficients in computing the shear stresses. This theory is used as a 
prime base in this finite element analysis. 

The finite element analysis of laminat ed composite plates has been presented by 
several authors, Reddy et. aL [Refe. 26,36], and Noor [Refe.37,38], for bending and 
vibration analyse. Mawenya [Ref 40], developed a general quadratic isoparametric 
muhilayer curved plate element, and Panda[Ref41], presented a superparametric, 
quadratic plate element for the plate bending analysis. Fortier and Rossettos [Ref 42], 
Sinha and Rath [Ref 43], analyzed firee vibration and buckling of thick plates of 
unsymmetric cross-ply construction. Dong [Ref 44] has given the solution for the 
dynamic response of a simply-supported rectangular plate. 

In the first part of the present work, a finite element model for a laminated conposite 
plate is developed using a simple, higher-order, shear deformation theory with 
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Hamilton’s principle for the formulation of the equations of motion [Refe. 45,46]. A 
standard, four node, rectangular element with seven degrees of freedom at each node is 
developed for the analysis of a flexible laminated plate having a constant thickness for any 
individual layer. The displacement model is so chosen because it can represent adequately 
the parabolic distribution of transverse shear stresses and the non-linearity of in-plane 
displacements across the thickness. A bending, free vibration and stress analysis problem 
is examined for different loadings, boundary conditions, and fiber orientation angles. The 
results are coiiq>ared with existing anafytical and numerical solutions. Hence the present 
element formulation demonstrates its applicability for a wide variety of laminate conposite 
plates. 

Given the result of continuous con:q)eting requirements for in:q>roving the weight, 
interdisciplinary performance, ten:q)erature stability, versatility, and reliability of propitlsion 
and aeroqjace conponents, the development of a new generation of composite materials, 
called ‘ intelligent/smart conq)osite materials ‘ is receiving growing attention, which is the 
concern of this dissertatioa 

The reader is referred to books, Cady [Ref 47], Dceda [Ref 48], and Nye [Ref 49], 
for piezoelectricity and for the development of the constitutive relations for piezoelectric 
materials. Tiersten [Ref 50], and Rogacheva [Ref 51], contain methods for solving the 
differential equations of the theory for piezoelectrical plates and shells, respectively. 

Several researchers have studied the interaction between the mechanical properties 
and the electric field. Crawley et. al. [Ref. 52], developed piezoelectric elements for 
laminated beams, and with his co-workers[Ref 53], has expanded the work and 
considered the piezoelectric actuators to be plies of laminat ed plate and used the Rayleigh 
Ritz method to study the deformations of a smart plate. They also modeled the induced 
strain actuation for a beam [Ref 54], and for a truss element[Ref 55]. Lee [Refe. 56, 
57], derived a theory for laminated piezoelectric plates based upon classical plate theory. 
His experimental [Ref. 58], results showed that PVDF or PVDF 2 actuators can generate 
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plate bending and twisting independently or simultaneously, and they are suitable for 
active d^Tnping control of a flexible structure. 

Tiersten [Ref. 50] modeled single-layer piezoelectric plates, including the charge 
equations, but did not study laminates. Wang and Rogers [Ref 59], applied the classical 
laminated plate theory to model laminated plates with spatially distributed actuators. 
Tzou et. aL [Ref 60], derived governing equations for piezoelectric shells using first order 
theory. Mitchell and Reddy [Ref 61], presented a hybrid theory for enhancing laminated 
plates based on modeling the electric potential through the laminate thickness with a 1-D 
finite element. Hagood et. al. [Ref 62], naodeled the effects of the dynamic coupling 
between a structure and an electrical network through the piezoelectric effect for a 
cantilever beam. State space models were developed for three important cases: direct 
voltage driven electrodes, direct charge driven electrodes, and an indirect drive case, 
where the piezoelectric electrodes were connected to an arbitrary electrical circuit with 
voltage and current sources. Ray et. aL [Ref. 63], presented the exact solutions 
for a two dimensional analysis of a plate composed of piezoelectric material. This study 
led to exact solutions for a con^xiske plate with piezoelectric actuators and sensors [Refe. 
64, 65]. Lee [Ref 67], developed an analytical approach via state space equations. 
Heyliger [Ref. 68], developed an exact solution for the free vibration analysis. 

The previous solutions do not provide the results for large complicated structures 
with integrated materials. Thus, the necessity for approximate techniques, such as the 
finite element method arises. A few papers have been developed addressing the analysis of 
intelligent structures by the finite element method. Allik and Hughes [Ref 69], presented 
a tetrahedral finit e element for three-dimensional electroelasticity. Based on this model, 
Tzou [Ref. 70], proposed a method for isotropic plates using isoparametric hexahedron 
solid elements. Kagawa et. al [Ref 71], developed a method for a transversely vibrating 
bar with electrostrictive transducers. McDearmon [Ref 72] and Tzou [Ref 73], 
developed a method for single cases of deformation. Ha et. al.[Ref 74] anal>^d a 
composite plate by using a three dimensional brick element. Both elements used in these 
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methods made the problem conqjlex, costly and required some special techniques to 
overcome inaccuracies and disadvantages of modeling a plate with 3-D elements. These 
elements display excessive shear stifBiess as the element thickness decreases. The two 
dimensional quadrilateral plate element developed by Hwang et. aL [Ref. 75], is more 
efficient than solid elements, but it appears to have restricted modeling capabilities. 
Hwang et. al. [Ref.76], proposed a model based on classical laminated plate theory, which 
neglects the effects of transverse shear stresses, and hence is inadequate for the analysis of 
moderately thick composite structures. Chandrashekhara, et. aL [Ref. 77], developed a 
model based on the first order shear deformation theory, which needed a shear correction 
coefficient. Static analysis or an intelligent plate was presented by Ray [Ref. 78], using a 
higher order shear deformation theory, which added additional dependent unknowns and 
made the problem costly to solve. 

In the second part of the present work, a finite element formulation is developed for 
modeling the static and dynamic response of laminated conqiosite plates with distributed 
piezoelectric actuators and sensors subjected to both mechanical and electrical loading. A 
sinople higher-order shear deformation theory with Hamilton’s principle is used to 
formulate the equation of motion of the system in which the piezoelectric layer is treated 
as a normal lamina. The displacement model represents the parabolic distribution of 
transverse shear stresses and the non-linearity of in-plane displacements across the 
thickness. The model is valid for both segmented and continuous piezoelectric elements, 
wWch can be either surfece bonded or embedded in the laminated plate. The piezoelectric 
layer can be isotropic or orthotropic, and the structure core can be anisotropic 
(graphite/epoxy, etc.) or isotropic (aluminum). A new, four-noded, bilinear, isoparametric, 
rectangular element with seven mechanical degrees of fireedom and one electrical degree 
of fi-eedom at each node is developed. The electric potential is treated as a generalized 
electric coordinate like the generalized displacement coordinates at the mid-plane of the 
actuator layers. The structural deformations due to electrical and mechanical loads are 
computed and compared to the available analytical and finite element resuhs. 
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The third part of this dissertation is directed toward shape control and optimization of 
the actuator voltages. Several works have been conpleted in vibration control. For 
exanple Baz et. al. [Ref 79], Challopadhyay [Ref 80], Chandrashekhara [Ref 77], Park 
eL at. [Ref 76], Tzou [Ref 81, 82], Anderson and Hagood [Ref 83], and Hagood et. al. 
[Ref 84]. Wang and Fuller [Ref 85, 86, 87, & 88], have done several papers for active 
structural acoustic control using piezoelectric actuators. A few studies have been done in 
the shape control and optimizatioiL Ghosh [Ref 89] showed a model for plate shape 
control by using PZT, Agrawal et. aL [Ref 90], developed a mathematical model for 
deflection using a finit e difference technique and estimated the optimal actuation voltages. 
Kirby [Ref 91], and Koconis et. al. [Ref 92], presented an anatytical method to determine 
the voltages needed to achieve a specified desired shape with minimum error between the 
actual shape and the desired shape. Based on Koconis method, Varadarajan et. al. 
[Ref 93], showed a model for shape control of a laminated plate. 

In the third part of the present woric, an optimization algorithm based on the finite 
element technique is developed to determine the optimal voltages applied to the actuator 
to minimize the error between the desired shape and the actual shape. The error flinction is 
defined as the mean square of the error between the points in the actual sur&ce and the 
points in the desired sur&ce for each element of the finite element grid. The original shape 
and the desired shape of the plate are specified. Thus the objective is to find the optimal 
actuator voltage to get minimum error between the desired shape and the actual shape. 
The analysis is based on small deformation theory, therefore, the specified desired shape 
must be within the region of small deformation fi’om the original shape. 

Two Matlab codes were developed. The first code, ‘COMPZ’ is an interactive finite 
element code. It is able to solve a laminated composite plate, with or without 
piezoelectric layers, subjected to mechanical and electric loads for different boundary 
conditions. The code is able to analyze either a conq)lete or a quarter plate analysis to 
save computational time, and it will perform a bending, fi'ee vibration analysis, and 
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a stress analysis for the individual layers. It is valid for different types of material for both 
piezoelectric or lamina ted layers ( anisotropic, isotropic, etc.). 

The second Matlab code, ‘OPTSHP’, is able to analyze the structure and confute the 
change in shape due to mechanical and electrical loads. This code includes a loop 
co ntaining an optimization algorithm coupled with a Matlab Optimization Toolbox to 
compute the voltage applied to each actuator to minimize the objective function. The 
numerical results are presented for each part of the analysis. 
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n. FINITE ELEMENT ANALYSIS OF A COMPOSITE PLATE 


A. INTRODUCTION 


The partial differential equation governing composite laminates of arbitrary geometry 
and boundary conditions can not be solved in closed form [Ref 95]. Analytical solutions 
of plate theories are available in the literature (Le., Navier solutions, Levy solutions). The 
Rayleigh-Ritz and Galerkin methods can also be used to determine approximate analytical 
solutions, but they too are limited to sinple geometries because of the difficulty in 
constructing approximation fmctions for complicated geometries. The use of numerical 
methods &cilitates the solution of these complicated problems. Among numerical 
methods, the finite element method is the most effective. 

There are several types of finite element models that have been developed for plate 
theories. These can be grouped into three major categories; displacement models, mixed 
and hybrid models, and equilibrium models. 

The displacement finite element models of plate theories are based on the principles of 
virtual displacements, where all governing equations are ejpressed in terms of the 
displacements. The mixed and hybrid finite element models are based on modified or 
mixed variations of the plate theories, in which both displacements and stresses are 
independently approximated. The equililsium models are based on the principle of virtual 
forces. Among the three types of models, the displacement finite element models are the 
most commonly used in commercial finite element progr ams . 

The objective of this chapter is to develop a finite element model for a laminated 
composite plate, using a single, higher-order, shear deformation theory. The model is 
then validated for different boundary conditions and loading cases. Bending, fi'ee vibration 
and stress analyses are performed using the proposed theory. 
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B. GENERAL THEROY 


Initial configuration 
of a vertical line 



Average deformation 
configuration of the 
vertical line 

Mid surface 

Actual deformation 
configuration of the 
vertical line 


Figure 2.1: Deformation of the normal to the mid-plane ‘Trom Ref [45].” 


A typical rectangular laminated plate has a length a, width b, and thickness t. The 
laminate is composed of a number of perfectly bonded orthotropic layers (laminae), which 
are placed sequentially one after another. A coordinate system is adopted such that the x-y 
plane coincides with the mid-plane of the plate, and the z axis is perpendicular to the 
plane. The present theory is formiilated based on the following displacement fields [Ref 
45]; 
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( 2 . 1 ) 


U{x,y,z) = U{x,y,Q) + z^,{x,ySi)-^z^Cxix,yS>) + z^g,{x>yS>) 

V(x,y,z) = Vix,y,0) + z^/x,y,0) + zV,(x, J'.O) + zV^(x,y,0) 
= Fo+z^^ + zV^ + zV^o 


iF(x,:v.z) = FF(x,>’,o) = iro 
where; 

C/q, Fgand are the displacements of a point on the reference sur&ce with coordinates 
(xo, yo, Zo). <f>^ and are the average rotations about the Y and X axes respectively of 

the normal to the mid-sur&ce of the undeformed plate, and z is the distance of a point 
from the mid-plane along the Z axis as shown in Figure (2.1). The remaining terms 
correspond to higher order rotations. Applying the conditions that the top and the 
bottom surfrice are free from transverse shear stresses; 

r„(x,>;,±r/2) = 0 and r^(x,>^,±r/2) = 0 (2.2) 

where t is the plate thickness. For an orthotropic plate, this means that the shear strains 
are zeros, Le. 


f„(x,y,±r/2)=0 and £^(x,>^,±r/2) = 0; (2.3) 

Differentiating equation (2.1) and substituting in equation (2.3) 
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= <^ / ^ + dfi^ / dc 

= + 3zV.(. / <*1.-1«! = 0 

sj = *>,. ± >Ct + (3 / + ^ / 4t = 0 

= <^ / <^+^ / 4^ 

^j®|z=±r/2 ~ ^ '*' ^ 4^l2=±//2 ~ ^ 

f>zlz=±r/2= ± + (3 ! 4)^V^ + ^l^ = 0 

By solving equations (2.4) and (2.5), we get 
Co = 0; Co = 0; 

Sio = -^^Uo + ^/^) and $-^0 = -^(<^ 3-0 + ^ / 4^) 
Substituting equation (2.6) into equation (2.1), we get 

U(x,y,z) = t/<, + z ^>,0 - (^,0 + ^ / ^) 

nx,y,z) = Fo + z Co - ^[7) (Co + ^ / 4 ^) 

lF(x,y,z) = W, 


(2.4) 


(2.5) 


( 2 . 6 ) 


(2.7) 
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The displacements at any point in the plate are given in terms of the seven unknown 
quantities A, and dt/v! 

C. STRAIN-DISPLACEMENT RELATIONS 


The linear strain-displacement relation in the Cartesian coordinate system is 




dUldc 





* = ^ 

{ajjc^^a^iac) > 



(aija+aviac) 



(3^fa+avi4>) 


By substituting equation (2.7) into equation (2.8), we get 




- 4(^wM^ + 



- 4(^w/ 4;^ + / 3r^ 


► ZZ * 

-4(2<?^w/<3:4; + ^,o.j- + /3^^) ' 

E„ 


4, - 4(7 / t)\a^ia + 4 )+ 



4 - 4(7 / r)^(^/4^+ 4 ) + ^ 0 / 4 ; 



(2.9) 
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where; 


= U ' 

= F • 

= Uq^ h- 9 

I 

£^ = Ao+f^oy> 

^yz ~ ^yO ■*■ 

Kl^Ky, 

f^y ~ ^y0^> 

■^jy ~ ^0,>’ ^y0,x’ 

Ki=^(^ldc + ^,,)lt^; 

=- 4 (^/ 4 ^+^yo)/^^; 
i:j=-4(^w/A^+^^,,)/(3f2); 

and 

4 = -4(2^w/at^+4<,^ + 4^)/3/“); (2.10) 

which includes linear strains^ curvatures, twists and other higher order curvatures. 

D. SRESS-STRAIN RELATIONS 

For anisotropic lamina, there are no planes of symmetry for the material properties. If 
there is one plane of material property symmetry, where the plane of symmetry is z=0 a 
material is termed monoclinic, and has 13 independent elastic constants . For a monoclinic 
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material, the stress-strain relationship is: 


o-/ 


■Qi 

Cn 

Q3 

Q4 

0 

0 ■ 

'X 

e. 



Qi 

Q2 

C23 

Q4 

0 

0 




C31 

Q2 

Qb 

Q4 

0 

0 

£2 

> 



C4, 

c« 

C43 

C44 

0 

0 




0 

0 

0 

0 

Q, 


^22 



0 

0 

0 

0 

Qs 

Qe. 



( 2 . 11 ) 


Since the normal stress Cz is small, it can be neglected. The corresponding 8z can be 
eliminated by putting Cz equal to zero in equation (2.11). This gives the reduced stress- 
strain relationship as 


O2 



C,2 

Q4 

0 

0 ‘ 

r 'x 

£2 



C21 

C22 

C24 

0 

0 

^y 


” = 

C4, 

C ,2 

Q4 

0 

0 


0 x 2 


0 

0 

0 

Qs 

^56 

^22 



0 

0 

0 

C« 

^ 66 . 

J 


where; 

Ml “Mi r ^ 
Ms 


Cj2 — Cj2 




23 


-^33 


Cj4 — 


Q3Q3 . 

c ’ 

M3 

c' ^ c 

Ml “Ml r ^ 
Ms 

r' ^ r ^^^Qs 
M2 ” M2 r" ’ 
Ms 


( 2 . 12 ) 
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^4 ' 

1 

II 

Q3Q3. 
Qa ’ 



^ 4 ,= 

= Q.- 

Q3Q3. 
c ’ 

^33 



C42 

1 

II 

C43Q3 

c ’ 

Ha 



c 

= <^ 44 - 

Q3Q3. 

■ C33 ’ 



C55 

II 


C5. 

il 

Qs 

= c • 

'-' 65 ’ 

and 

4 

0* 

II 


(2.13) 


or in short form; 

0/ “ ^3^73 ^ Q 3 

for 

i, j = 1,2,4 


II 

for 

ij = 5,6 


and the QcoeflScients [Refe. 96,97]; 


Qi ~ 




■"13 


1 

E 2 E 2 A 


C 23 - 


C33 - 


E2E2A 

£,£3/1 ’ 

t 23 ,+U 2 ,L >32 

^3 + ^ 2 l ^23 

E2E2A 

E^E 2 A ’ 

l-u, 3 t; 3 i 

£,£3^ ’ 


L >32 + L>, 2 U 3 , 

1223 + t* 2 ,Ui 3 

E^E^A 

£,£2^ ’ 

\-O^^V 2 ^ 


£,£2^ ’ 



(2.14) 
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and; 


Qs ~ ^13» 

Qe “ ^23» 

1 - - ^ 1^3 - 

£,£ 2^3 

wiiere £, ,£2 ,£3 are the Young’s moduli in the 1,2, and 3 directions, respectively; is 
the Poisson’s ratio for transverse strain in the j-direction >^n stressed in the i-directioa 

For orthotropic lamina there »e two orthogonal planes of material property symmetry 
for the material, and symmetry will exist relative to the third mutually orthogonal plane. 
The stress-strain relations equation (2.11), in coordinates aligned with principal material 
directions have no interaction between normal stresses and shearing strains; Le., 
C ,4 = C 24 = C 34 = C 4 J = 0, Thus the stress-strain relation of an orthotropic lamina is; 

fc;, c;2 o o olf^/ 

<^y C^l C 22 0 0 0 Sy 

0 0 C 44 0 0 (2.15) 

O’* 0 0 0 Cj5 0 

[0 0 0 0 

If the fibers are oriented at an angle 0 with the x-axis as shown in Figure 2.2, the 
transformed stress-strain relation for a lamina will be 
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Figure 2.2: An off-axis unidirectional lamina. 


r ^ 

O’, 


■Qu 

Qi2 

Q >4 

0 

0 ■ 


f ^ 



Q 2 , 

Q22 

Q24 

0 

0 




► — 

Q 4 . 

Q42 

Q44 

0 

0 



o-» 


0 

0 

0 

Q55 

Q56 





0 

0 

0 

Q 65 

0 . 66 . 




in short form; 

W] = \q\{£] 

where; 

Qii = Cjicos^ 6 + 2(Cj 2 + 2C44)cos^^sin^ 6 + Cjj sin^ 6 
Qi 2 = (Cii + ^22 - AC^J^)coi6 sin^^+ Cj2(cos'‘^-i- sin^ 6) 


(2.16) 


(2.17) 


(2.18) 



Q,4 = (C^i - 2C44 - Cj2)cos^0sin^+ (Cjj - + 2C44)cos^sin^^ 

Q22 = <^isin^^+ 2(C^2 +2C44)cos^ftm^^ + C^cos^^ 

Q24 = (^1 “ 2C44 - Ci2)cos^sin^^+ (0^2 - Cjj + 2C^^)oos^6&m6 
Q44 = (Cii + “ 2C^44 )cos^ 044(00$'’ + sin^^) 

Qj 5 = Cjj cos^ 6 + C^sin^^ 

Q54 = (Cjj - C^)cosftin^ 

Qj4 = Cj5sm^^+ 


E. STRESS RESULTANT-STRAIN RELATIONS 


Combining equation (2.9) with equation (2.16) and integrating layer-by-layer over 
the thickness, the following relations of the stress resultant are obtained [Ref 96]; 


= £%,(l,z,2')dr; (2.19) 

iNy,My,Py)= [^^^tTy(l,Z,Z^)ct; 


The above equations can be set in the matrix form as; 
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or in contact form 


(W) =[!)]{«•) 


^ere; 

[n]=Rigidity matrix 

[a] =Extensional stiflfiiess matrix 

[ 5 ] = Coupling stif&iess matrix 

[D] = Bending stif&iess matrix 

[£l,[F], and[H] = Higher order matrices 

which are given for n layer by: 


t f" Q,(ll,2‘,2>yy)d2 

*=1 * 


(forj,i=l,...,6) 


( 2 . 20 ) 


( 2 . 21 ) 
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F. VARIATIONAL PRINCIPLE 


The generalizd fonn of Hamilton’s principle, vsiiich states that the variation of the 
Lagrangian during any time interval must be zero, has the form: 

J[^^<J(3 + irVr = 0, (2.22) 

where and t^ are the time interval. 

The Lagrangian 3 of a body is defined by the summation of all kinetic energy X and 


potential energy ti . 

3= [{%-K)dV 

(2.23) 

where; 


(2.24) 

* = 

(2.25) 


Thus the Lagrangian is 

3 = dV (2.26) 

\^ere; 

^ is the velocity (time derivative of the displacement q); 

3 is the Lagrangian; 

V is the volume. 

The virtual work SfV done by the external applied force is 

SfV ^ [[Sq]'{P,)ds (2.27) 


where; 


s is the surfece areas at which the mechanical load is applied 
P, is a surfece load vector (N/m^). 

Substituting equations (2.26) and (2.27) into the Hamilton’s equation (2.22), yields the 
variational equation as: 


Since all the variations must vanish at r = t, and / = /j, by substituting equation (2.17), 
and taking the variation, the variational equation takes the form: 

lp{SgY{g}<iV*l{&Y[Qle)dV-l{s,Y{P,}ds = 0 (2.29) 

G. FINITE ELEMENT FORMULATION 

The objective is to define the degrees of fiwdom Uo,Vo,<t>xo,<|>yo,Wo,Wx, and Wy in the 
plate in terms of nodal displacements and rotations by using a bilinear, isoparametric, 
rec tang ular element with four node. Each node of the element has seven degrees of 
fi’eedom. S imilar interpolation functions for the element coordinates x and y; the in-plane 
displacements Uo and Vo and the two rotations <|>xo and (((yo were used. They were defined 

by; 

p = 

i=l 

where; 

p is the value of the variable at any point in the element 
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(Xc-a/2,yc-b/2) 

(Xc-a/2,yo+b/2) 

(-1,*1) 

(-14) 

Actual element 


Master element 



^ = 2(x-x^)/fl, Tj = 2{y-y^)lb 
Figure 2.3: Actual and master element. 

Pi is its value at node point 

Ni is the interpolation function, which in the natural coordinate system is; 

Ni = 1 / 4(1 + 44 ^X 1 +m) (2.31) 

where; 

^ and Tj are the local coordinates of the point, and 4 = -1,14-1 and 
^ = -1,-144 for / = 1,..,4, as shown in Figure (2.3) 

The transverse displacement is interpolated using a non-conforming shape function [Ref. 
98], which can be given by; 
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w(x,:v) 


[/, ft \ A g2 K A ft ^ /a ft 



(2.32) 


where for / = 1 -► 4 

/ = (1 /8) (1 + + Vn (2-33) 

g, = (a/16)^(l + -1) 

= (h/16);z (1 + mf<y + ^4X/7?7, -1) 

where; 

^ = 2{x-x,)l a, t}-2(j- y^)/b, (2.34) 

are the centroid coordinates of the plate, and the non-dimensional coordinates of 
the nodes are (-1,-1), (1,-1), (1,1), and (-1,1), respectively. 

The nodal displacement vector at the first node point on the reference surfece is; 


{ft} “ [^01 IqI ^01 ^.JTl ^J'l] 

(2.35) 

and the element displacement vector \q^ is; 


{ft} = [ft ft ft ftf 

(2.36) 

The generalized displacement vector at a point is 


{9)=[^0 K 4>yO % Wx ^.y] 

(2.37) 


Substituting equations (2.30) and (2.37) into equation (2.10), we obtain the strain vector 
at the mid-plane in a matrix form 
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The first, second, and higher order derivatives of the interpolation functions \^ch are 
expressed in equations (2.38) and (2.39) can be evaluated with respect to the global 
coordinates. An add i tional conqjutations were involved ■winch can be eiqiressed as. Let ^ 
represent one of the interpolation fiinctions which has been used (Le. , _/], gj, and/or 

h)- 

The first order derivatives with respect to the global coordinates are related to those with 
respect to the focal (or element) coordinates according to: 






-1 














A 














(2.40) 


where the Jacobian matri x [ j] is evaluated using the ^proximation of the geometry wiiich 
may take the form: 


or 


x=x, + ^a/2; 

y = yc + i7b/2 


x = 

y = 


X1+X4 

2 

yl•^y4 

2 


+ ^a/2; 
+ i7b/2 


(2.41) 


The second order derivatives of Ti with respect to the global coordinates (x, y) are 
given by: 
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(2.42) 




f 

\^%] 




-N" 

i 

a^% 

af 

-[j=} 

. ^ . 

A^. 


\ 





where; 








\dt]) \^) 

dj 

dqd^. dr]d^ 




2 — 
dr\ dri 

drj dt] 


w= 





ae 

^x 

£y 



^x 


dnd^ 

a^an 


(2.43) 


(2.44) 


Tte elements of the matrices [j,] and [jj] 


are computed using equation (2.41). 


The generalized strain at a point is related to the reference surfece str ain as: 


where; 


(2.45) 
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(2.46) 



lOOOOzOO 
0 1 0 0 0 0 z 0 
0 0 1 0 0 0 0 z 
0 0 0 1 0 0 0 0 
0 0 0 0 1 0 0 0 


0 

0 

0 

z^ 

0 


0 z^ 0 

0 0 z^ 

0 0 0 

0 0 0 

z^ 0 0 


0 

0 

z^ 

0 

0 


The reference surfece (mid plane) strain can be expressed as; 

Thus; 

By using the transformed stress-strain relation, 

io’]=[0]{4 

in equation (2.49), 

W]-lQlBpU} 


(2.47) 


(2.48) 


(2.17) 


(2.49) 


H. EQUATION OF MOTION 


By substituting equations (2.48) and (2.49) in the Hamliton’s equation (2.29) gives 
the ^^em equations of motion: 


l{6q.}’ANY[NW{i.]* l{Sq^’[BY[B,"(\Q\B\B\dv{q,]- Vl'fn}* = 0 

(2.50) 
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In matrix form, the equatioDS of motion can be e^qyressed as: 


+[<;]{,,) = {/v) 

>^i)ere the element mass matrix is; 

\K]= il{NY\NW 


The sh£q)e function matrix [AT] is given by [Ref 47]; 





0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

o' 


0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 


& 

h 

0 

0 

0 



K 

0 

0 

0 





I 


denotes the matrix assemblage for the element four nodes, and 


(2.51) 


(2.52) 


(2.53) 


N = 


/. 0 
0 /, 
0 0 
0 0 
0 0 
0 0 
0 0 


0 

0 

h 

0 

0 

0 

0 


0 

0 

0 

h 

0 

0 

0 


0 0 0 

0 0 0 

0 0 0 

0 0 0 

/, 0 0 

0 /j 0 

0 0 /j 


(2.54) 
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in /j and I 2 flre the normal and rotational inertia, respectively. 


/.I ' 


(2.55) 


\niiere p‘ is the mass density of the i* layer. The element stif&iess matrix has the form; 


= llBflDlB\dA 


(2.56) 


where the rigidity matrix [d] for n number oflayers is given ly; 


[a] = Z Q^{U,^. 2 \z\z^)dz (for j, i =1,.. .6) 


(2.57) 


The rigidity matrix D is given in equation (2.20), and the mechanical excitation force is 
given by; 


Thus the consistent load vector at node i is 

O' 


“ I ^0 


0 

0 

0 

8i 

L*(, 


Ids 


(2.59) 


where; is the intensity of the load per unit area. 
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Assembling all elemental equations gives the global dynamic system equation: 


[w){«}+[a:]{?} = {/>„) (2.60) 

I. NUMERICAL INTEGRATION 

The elements of the equation of motion (2.S1), which are giv»i individually in 
equations (2.52), (2.56) and (2.59) can be integrated numerically using Gauss 
Legendre quadrature (see App. A). A M integration technique of 3x3 Gauss points 
is used to perform the integrations, and it gives satis&ctory results (comparing to the 
exact and available finite element solutions) for all plates. 




[*.]= 


i«i 


0 

0 

0 

0 

f> 


gf 


\J\d^Tj 


(2.61) 


(2.62) 


(2.63) 
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J. FREE VIBRATION ANALYSIS 


A free vibration analysis can be performed using the equation of motion (2.60), 
where the right hand side equals zero (the work done by external forces is equal to zero 
for the free vibration analysis). The effects of the degree of orthotropy, of different 
stacking sequences, aspect ratio (span-to-thickness ratio) of the structure, numbers of 
layers and orientation angles of the fibers, on the natural frequencies and mode shapes 
are demonstrated in the validation part. 

K. STRESS ANALYSIS 


Snoonwl sK«» ftfdisn 



Figure 2.4; Local discrete least squares smoothing “From Ref. [99].” 

By solving equation (2.60) for a static deflection, the element displacement vector is 
computed , v^diich can be substituted into the generalized strain at a point equation 
(2.48), then into equation (2.17) to get the stress vector at any point in the structure. 
Stresses are obtained at 2 x 2 Gauss sampling points, as shown in Figure (2.4) [Ref 99]. 
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In local discrete smoothing [Ref. 99], it is assumed that is an exact least 

squares fit to selected values of at the Gaussian integration points sho^vn in 

Figure (2.4). (exact least squares fit is one ^\1uch passes through all the sanpling points 
and is therefore interpolatory in nature). Thus, the smoothed comer nodal stresses 
CT, ,CTj, and a^ may be obtained fium the ejquession. 



1 ^ 

1 + 

1 


1 


2 

“2 

2 

2 


1 


1 



~2 

2 

"2 

2 

^// 


1 


1 


2 

~2 

2 

2 

^IV 

2 

>4 

_1 

~2 

1 

T . 



(2.64) 


where and are the unsmoothed stresses at the Gauss points 

(cTj = as shown in Figure (2.4). The smoothed stress values are then modified 

by finding the average of the nodal stresses of all elements meeting at a common node. 


L. VALroATION: 


To test the validity of the finite element analysis technique and to establish its range of 
£q>plicability, numerical examples are investigated. 

Exanqjle 1. 

A three- layer (0°-90®-0°), simply supported, square plate of three different span to 
thickness ratios (X.=10,20,100) is tested . The total thickness of the 0° and 90® layers is 
the same. Layers at the same orientation have equal thicknesses (i.e. ti=t 3 and ti+t 3 =t 2 ). 
The material constants are as follows; 
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Table 2.1: Boundary Conditions for various different conditions (x - constant). 


Sin^fy supported 

Clanped 

Free 

U^O; 

F = 0; 

(7 = 0; 

F = 0; 

C/9t0; 

V^O; 


*, = 0; 

^<0 ~ 

II 

0 



1^ = 0; 

* 0; 

w = 0; 

35 
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p 

_ ! 

w^O; 

w^9t0; 


w^ = 0 



En = 172.4 GPa (25x10* psi); 

E 22 = 6.9GPa(10*psO; 

Gi 2 = Gi3 =3.45 GPa (0.5x10* psi); 

G 23 = 1.38 GPa (0.2 X 10* psi); 

V 12 = Vi3=0.25 

Aspect Ratio( length to thickness ratio)=10,20 and 100; 

The plate is subjected to a double sinusoidal load 
q=qosin(7tx/a)sin(jty/b) 

where a*b=L; the plate lengtL Table 2.1 gives the boundary conditions on the edge 
X = constant. S imilar statements can be made for the edge ^ere y is a constant by 
interchanging the subscripts for x and y. A 4x4 mesh quarter plate system with 3x3 
sampling points, is used in all exaiEq)les. The results are given in Table 2.2 in normalized 
quantities. The results are normalized as follows: 

= (1 / ; 

= (1 / 

w = n^Qyi / 12X* tq ^; 
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Table 2.2: Three-Layers cross ply (0°/90®/0®) square plate [Ref. 45] subjected to sinusoidal 
loading (r,=r 3 = r/4, fj" 
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wiiere; 

S is the reference source A Present finite element. 

BFE2/PandaandNatarajan[Ref. 41]. C FEl/Reddy [Ref. 27] 

D FE3/Mawenya and Davies [Ref 40]. E ISPQ/Moser et. aL [Ref 104]. 

F FE4/Phan and Red<fy [Ref 33]. G Exact/Pagano and Hatfield [Ref 16]. 

CPT Classical plate theory. 

and 

X = alt, vsiiere r is the structure thickness. 

Q = 4Gj2 + [£„ + £ 22(1 + 2 ^ 12 )]! 0 - ^2t>2i); 
yv^=w{al2,al2,0y, 

CTj = ^(fl/2,fl/2^1/2) 

cTj = o>(a/2,fl/2^1/4); 

<73 = (0,04:1/2); 

- error in w fi-om exact value 
= error in tr, fi’om exact value. 

E„ = error in a, firom exact value. 

E„ - error in a, fi’om exact value. 

Table 2.2 shows that the numerical values of all quantities are converge to the CPT 
results by increasing the span to thickness ratio k . A model devebped by Reddy [Ref 
27], based on the penahyAfNS, (YNS is a theory of Yang, Norris and Stavsky) with eight 
node quadrilateral elements and five degrees of fi«edom per node is perform better for 
small span-to-thickness ratios k . Ebments developed by Panda, Mawenya and Moser, 
either have too many structural degrees of freedom, or have unsatisfectory performance 
conqjarison to the present model. The Phan and Reddy element gives very good results 
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conqjared to the exact solution for the thick plate (small X), but the error of this element 
increases with the increase the span to thickness ratio. Thus the conclusion is that for 
higher span to thickness ratios one can use the present element, and for lower span to 
thickness ratios, Phan’s element can be used for better accuracy. For length to thickness 
ratios equal to fifty. Figures (2,5), (2.6), and (2.7) show the different stress distributions 
for a composite laminate using the present finite element method. 
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Normal Stress distribution Of Simply Supported Composite Plate 



Figure 2.5: Distribution of the in-plane normal stress Oxx for a singly supported 
laminat ed square plate (0*’/90V0°) subjected to double sinusoidal load. 
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Normal Stress distribution Of Simply Supported Composite Plate 



Figure 2.6: Distribution of the in-plane normal stress Cyy for a singly supported 
laminated square plate (OVmVO®) subjected to double sinusoidal load 
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Shear Stress distribution Of Simply Supported Composite Plate 



Figure 2.7: Distribution of the transverse ^ear stress Oxz for a singly supported 
laminat ed square plate (0°/90®/0°) subjected to double sinusoidal load. 
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Exanqjle 2. 


The fallo wing mater ial constants are used in the free vibration analysis of a composite 
plate [Ref. 46]; 


E,/E2 = 40; Gi2=Gi3 = 0.6E2; 

G 23 =0.5E2, vi2 = 0.25 

The dime nsionless fundamental frequency of the present finite element is con:^)ared 
with different methods of analysis given in Table 2.3, for span to thickness ratio equal to 
five. The effect of span to thickness ratio and the fiber orientation angles of the 
din^nsionless fundamental frequency are given in Figure 2.8. These results validated the 
finite element code. 

Based on the results, the fundamental frequency can be increased by increasing 
lamination angle up to 45° except the case of two layers in which it decreases. Increasing 
the number of layers with fixed thickness increases the fundamental fi^uency due to the 
decreasing of flexural extensional coupling effect. The dimensionless frequency increases 
as long as the span to thickness ratio decreases and reaches a maximum value at a fiber 
angle equal to 45°. 
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Table 2.3: Dimensionless fundamental frequency for four layers (0°/90®/0°/90®) with span 
to thickness ration equal five 


Method 

Normalized Natural Frequencies’ 

Present FE 

0.4500 

HSDT* 

0.44694 

HSDT" 

0.44686 

HSDT' 

0.44686 

FSDT" 

0.44083 

FSDT' 

0.45083 

Noor" 

0.42719 


0.66690 

CPT' 

0.66690 


Where; 

HSDT is the higher order shear deformation theory [Ref 35]. 

FSDT is the first order shear deformation theory [Ref 35]. 

CPT is the classical laminate plate theory. 

* Results obtained using the finite element solution [Ref 35]. 

Results obtained using the Navier Solution. 

‘ Results obtained with the exact solmion [Ref 35]. 

** Results obtained by applying a finite difference scheme to the equations of the three- 
dimensional elasticity theory. 

* “FromRef [35]” 


44 











Wn=w*sqrt(row*t''2/E2) 


Effects Of Fiber Angles & Span to Thickness Ratios On Dimensionless Natural Frequency 



Figure 2.8: Effects of the fiber angles & span to thickness ratio on the dimensionless 
fundamental natural frequency for a laminate (0°/-6®/0°/-0°). 









m. FINITE ELEMENT ANALYSIS OF A SMART 
COMPOSITE PLATE 


A. INTRODUCTION 

The rapid development of high-speed con 5 )uters has fecilitated the use of 
computational techniques in a variety of engineering applications. The finite element 
method is one of the most popular and powerful techniques in modem engineering design 
and analysis of con^licated structures and multifield problems. Most of the research on 
active piezoelectric stractures has focused primarily on ejq)erimental and theoretical 
studies, and there has been little development of general purpose piezoelectric finite 
element codes [Ref 99]. 

In general, e;q>erimental models are limited by sizes, cost, and other laboratory 
unknowns. Theoretical models can be more general, however, analytical solutions are 
restricted to relatively sin 5 )le geometries and boundary conditions. In the case of 
complicated geometries and/or boxmdary conditions, both theoretical and e)q)erimental 
techniques encounter technical difficulties. Thus, finite element development becomes 
very inqjortant in the modeling and analysis of a elastic/piezoelectric coupled sy^em. A 
typical elastic/piezoelectric stmcture is con^sed of a main elastic continumn, such as 
aluminum or graphite epojty, with coupled (surfece bounded), or embedded piezoelectric 
sensors and actuators. The thickness of the main stmcture can be about two or three 
times thicker than that of the piezoelectric layer. Thus, it would be very inefficient and 
time-consuming if the entire stmcture were modeled by isoparametric hexahedron or 
tetrahedral solid elements. Thus, the development of a new finite element model applied 
to piezoelectric con:qx)site laminated plates, based on a simple, higher-order, shear 
deformation theory is essential to investigate the stmctural response due to the applied 
field and ways to control the shape of the stmcture. 
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In this cb^ter a finite element model is developed and a single higher-order shear 
deformation theory derived with Hamilton’s principle used to formulate the equations of 
motion of the structure. The model is valid for a piezoelectric layer either surfece bonded 
or embedded in the laminated, plate. The piezo-lamina could be patches or completely 
cover the sur&ce as shown in Figures (3.1) and (3.2). A four noded, bilinear, 
isoparametric rec tangular element with seven mechanical degrees of fireedom and one 
electrical degree of fi'eedom is developed. The electric potential is treated as a generalized 
electric coordinates, similar to a generalized displacement coordinates at the mid-plane of 
the actuator layer. 
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b. Simply supported plate with piezoelectric layer cmnplete cover the surfece 

Figure 3.1 Composite plate with piezoelectric sensor and 
actuator conpletely covering the surface 
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Piezoelectric Actuators 


Piezoelectric Sensors 


Figure 3.2: Con^site plate with distributed segments piezoelectric 
actuators and sensors 
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B. VAWATIONAL PRINCIPLE 


The linear piezoelectric constitutive equations coupling the elastic field and the electric 
field can be e3q)ressed as [Ref. 47]; 
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(3.1) 


(3.2) 


Equation (3.1) describes the direct piezoelectric effect, which means a charge/vohage 
generated by an imposed force/pressure to a piezoelectric material. The converse 
piezoelectric effect is described by equation (3.2) in which induced stress/strain are 
induced due to an externally applied vohage/charge. The equations can be written simply 
as: 


{D}=H^{4+MU} (3.3) 

M=[c]{s}-[e]{E} (3.4) 
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where; 

{ D) = Electric displacement vector (Coulomb / meter square) 

[e] = Dielectric permittivity matrix (Coulomb / meter square) 

(f) = Strain vector 

[f"] = Dielectric matrix at constant mechanical strain (Farad / meter) 
(permittivity conqwnent) 

[ E] = Electric field vector (Volt / meter) 

{<t} = Stress vector (Newton/ meter square) 

{c} = Elasticity matrix for a constant electric field (Newton / meter square) 


It is assumed that the principal material coordinates coincide with the coordinates of 
the problem being analyzed. Thus the constitutive relations for a material having 
orthohombic mm2 symmetry, including piezoelectric effects, are given by: 
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In this analysis, the plane stress approximation is made by setting <7^ = 0, the strain 
is eliminated from the constitutive relations, which will takes the following form: 


a' 
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where; 

- _ _ . 

®31 “ ®31 „ ®33> 

*'33 

— ^23 

®32 ~ ®32 ” ®33» 

*'33 

^24 ~ ^ 24 ’ 

^15 ~ ^15» 

^22 ~ ^2’ 

and 


^3 “ % ^ 

<^33 


(3.7) 


(3.8) 
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Table 3.1 Analogy between mechanical and electrical quantities. 


Mechanical 

Electrical 

Force density Ps (vector) 

Charge density p (scalar) 

Displacement q (vector) 

Potential O (scalar) 

Stress c (second -order tensor) 

Flux density D (vector) 

Strains (second-ordertensor) 

Electric field E (vector) 


The elastic coefficients Cjj and Cjj can be obtained in the same way as Cjj and Cjj, which 

are given inequations (2.13) and (2.14), and the transformation (ifnecessary) can be 
confuted by the same method described in equation (2.18). The analogy between the 
mechanical and electrical quantities [Ref. 69] is given in Table 3.1. 

The Lagrangian, 3, of a piezoelectric body is defined by the summation of all kinetic 
energy ,X, and potential energy ,ti, (including strain and electrical energies): 


:i=[{X-n)dV (3.9) 

where; 

= (3.10) 

and 

» = (3.11) 

Thus the Lagrangian is 




(3.12) 
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w^ere; 


q is the velocity (time derivative of the displacement q); 

3 is the Lagrangian; 

V is the piezoelectric volume. 

The virtual work, SfV, done by the external surfece force and the applied surface 
charge density, //, (CW) applied to the piezoelectric body is; 


(313) 

where; 

5, and $2 are the sur&ce areas at which the mechanical and the electrical loads are 
applied; respective^ 

P, is a sur&ce load vector (N/m^) 

<P is the electric potential (volt) 

The minus sign in equation (3.13) occurs because in the variational principle for the 
electromechanical medium it turns out that the electric enthalpy H = U- ; takes the 

place of the internal energy function i/ in the Lagrange density; Le., the effective 
electrical energy content of // is opposite in sign to that of U . 

By using Hamikon’s principle 

l's(:i+W)dt = 0, (3.14) 

where /, to tj is the time interval, and all variations must vanish at t = r, and t = tj. 
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Substituting equations (3.12) and (3.13) into equation (3.14) yields the variational 
equation as: 




By substituting equation (3.3) and (3.4) and taking the variation gives: 


(3.15) 


i: 


lU) ++ {^]y]{E]w 


Wr = 0 


(3.16) 


Since all the variations must vanish at r = and t = t 2 , the variational equation takes the 
form: 


i\A^Y{i} * - (&rw{£} - {«r^r{«} - {m'U']{E]yv 
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C. FINITE ELEMENT FORMULATION 


To generate the electro-elastic matrix relations for a finite element, it was assumed 
that the surfiice of the piezoelectric layers which are in contact with the laminated 
substructures are suitably grounded. Also, since the thickness of the piezoelectric layers is 
very small, it is reasonable to assume that the electric potential fimctions, \s1iich yield zero 
potential at the interfece between the actuator and/or the laminated substructure and 
provides a linear variation across the thickness of the sensor or actuator layer are as 
follows: 

<t>\x,y,z)^{z-h^)<S>l{x,y) (3.18) 

where can be treated as the generalized electric coordinate similar to the generalized 
displacement coordinates at the mid-plane of the actuator and sensor layers. The 

generalized electric coordinates at any point within the element can then be e^ressed in 
terms of i nodal variables values via interpolation fimction N^: 

= [A'.K®;) (3.19) 

where is the nodal generalized electric coordinate vector and is given by 

{<«>.'}=[<!>.• K K (3-20) 

with 0QI (i = 1,..,4) is the generalized electric coordinate at the i* node of the element 

and is the shape fimction matrix, that is, the same linear shape fimction used for 
element coordinates x and y, which are given in equation (2.31). 

iV,= 1/4(1+ #4)(1 + 7 ;x) (3.21) 

where, 

= N, N, W,] (3.22) 
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Thus; 




(3.23) 


The electric field vector {£} is defined by the electrical potential energy d> by using 
gradient operator V: 

{£}=-Vd>^ (3.24) 


Substituting equation (3.23) in equation (3.24), gives the electric field vector; 


where; 


and 


[E] = 


[s.] = v[w.] 
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_ 1 



N, 



0 0 

-iz-h^) 0 

0 -1 


(3.25) 

(3.26) 


(3.27) 


The strain vector {f®} at the reference surfece is defined by; 

{^"}=[5]{9,) (3-28) 

where; 

\q^ is the element displacement vector given in equation (2.37). 
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[ B] is the nodal strain-displacement matrix given in equation (2.38). 

The generalized strain at a point is related to the reference surfece strain as; 


{^)=[5,F! (3,29) 

Substitute equation (3.28) into equation (3.29) yields 

(3-30) 

where; 

lOOOOzOOO Oz^O o' 

OlOOOOzOO 0 Oz^O 
0 0 1 0 0 0 0 z 0 0 0 0 z^ (3.31) 

00010000 z^ 0 0 0 0 

000010000 z^ 00 0 

Substitution of equations (3.25), (3.30) into equation (3.17) yields; 



I {«,)>]"{ P,}*, + 1 {»;) V.f (r - =0 


(3.32) 



D. EQUATIONS OF MOTION 


1. Plate With Actuators Only 
The variational equation (3.32) can be written in the form; 

[AN\\N\dV%]{SiX * 

(3.33) 

-1 [Af]'{?,}*,{<%,)'+ {«>'.)' = 0 

Thus the equation of motion in matrix form is 

+[^w]ke} ~ (3-34) 

[*,.,]{«,)+M{4-:)={«) 

After adding artificial, linear, viscous damping to equation (3.34), the equations of motion 
will take the fonm 

[M.]fe} + K,]fe} + [^,,]{9.}-[^^]{‘I>o} = {Pm} (3.35) 
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where; 


\m.\ = [AnXWv 

“ ljiNflm1iN]d(area) 


The shape function matrix [A^] is given by 





0 

Nj 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Nj 

0 

0 

0 

0 

0 

Nj 

0 

0 

0 

0 

0 

y; 

& 

A, 

0 

0 


8i^ 


0 

0 

fi^ 

Sijf 



and 


[m] = 


0 

0 

0 

0 

0 

0 


0 0 0 

/, 0 0 

0 0 
0 0 /j 

0 0 0 

0 0 0 

0 0 0 


0 0 0 

0 0 0 

0 0 0 

0 0 0 

/, 0 0 

0 /j 0 

0 0 /j 


in which /, and ^ are the normal and rotational inertia, respectively; 




where p' is the mass density of the /'* piezoelectric layer. 


(3.36) 


(3.37) 


(3.38) 


(3.39) 
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The Hamping matrix |c^ J is defined as a proportional danping, i.e. 


and 

(3.40) 

[*«]-iisnArnAj«]* 

where the rigidity matrix [d] for n number of layers is given by 

(3.41) 

[^j] = S f ' Cj,(l,z,z^z^ 2 ^ 2 ‘)a!z (for j, i =1,...6) 

and its elements are given explicitly in equation (2.20) and 

(3.42) 

[*^] = [[BY[B^[t\z\B^\iV 
= l\BY[Be2\B,)iA 

where for m number of piezoelectric layers is given by; 

(3.43) 

*=i ^ ■' 

and 

(3.44) 

[*-]==[*^]' 

(3.45) 
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and 
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To specify the voltage, the distributed applied charge density required to create a 
prescribed electric potential distribution on the surfece of the actuator layer can be 
determined as; 


= - 


4 








(3.52) 


where />c^and ^y^are the distances from structure middle surfece to the outer and the inner 


sur&ces of the piezoelectic actuator layer, respectively. Assembling all the equations 
gives the global dynamic system equations 


[M]{9) +lC]{^i +[JC„]{«) -[i,.]!*} = {F) (3.53) 

Note that the mechanical equation is coiq)led with the electrical equation; in which {F} is 
the global mechanical excitation and {G} is the electrical excitation. In active vibration 
control application, {G} is the feedback voltage determined by the control algorithm. 
In the static case, equation (3.53) becomes; 

By performing a condensation of the {d>} degrees of freedom from equation (3.54), the 
equation of motion will be as follows: 

lr]{4) = {F} (3.55) 

where; 

[r] = [x„] + [vI«~r'K] (3.56) 
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and the electrical potential vector is; 




(3.57) 


v^diich can be used to calculate the voltage distribution. In the free vibration analysis, { G} 
is set to zero so that the voltage distribution associated with each mode can be estimated. 

Also from equation (3.55), the feedback force jiy | can be expressed as 

P.58) 


2 . Plate With Actuators And Sensors 

For actuator and sensor the variational equation (3.32) will takes the form: 

lAN'([N]dv{q,]{SqX 

- I [A'f {?,}&,= 0 
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(3.59) 



Since dq^ and are arbitrary, equation (3.61) is satisfied only if 


= {/■„} ( 3 - 60 ) 

[t«,Ik} + [*Ml{®;L={«! ‘Actuator’ 

[*•»!{«.} +[*«]{<*'o), ' » 

After aHHing artificial damping and assembling all the equations, the global dynamic 
system equations are; 

[M]{q] +[C]{^} +[/:„]{?)= {F} (3.61) 

[A:c.,]j^}+[A^L{o}„ = {G} ‘Actuator’ 

+[^l{^}a = 0 ‘Sensor’ 

where {g}, {0}„, and {<!>}, are the global, nodal, generalized, displacement coordinates, 
and the global, nodal, generalized, electric coordinates for actuator and sensor, 
respectively. 


66 



The global, nodal, generalized, electric coordinates can be condensed using equation 
(3.61), yielding the system equations. 

(i/]{«j+[Cl{j)+[r]{,) = {f'} (3.62) 

^^ilere; 

{r} = {F) + [jc^] (3.63) 

Since the charge applied to the sensor layer is zero ( G=0), the voltage from the sensor 
layer can be written as; 

( 3 « 5 ) 

The values of the matrices in equation (3.60), m^, , k ^^, and can be 

conqjuted from equations (3.36), (3.41), (3.43), (3.45) and (3.46). The matrices , 
and k^^^ can be computed the equations (3.43), (3.45), and (3.46) for the sensor 
layers. The mechanical and electric loads are confuted using equations (3.49) and (3.51). 
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E. NUMERICAL INTEGRATION 


By using Gauss Legendre quadrature, a full integration technique of a 3x3 Gauss 
points (see Appendix A) are used to con:q)ute the integration of each element of the 
equations of motion of both cases, using actuators only or using actuators and sensors 
(3.34) and (3.60). Thus equations (3.36), (3.41), (3.43), (3.45), and (3.46), can be 


nxunerically integrated as follows: 

l[NX[mlN\iA. (3.66) 

(3.67) 

[t^] = I [BY[BeZ\B^]u, (3.68) 

[i«,] = £ [i!»]'[Z^'Z][s.]<i4, (3.69) 
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0 
0 
0 
0 

f, 

gi 

h 

'o' 

0 
0 

W't f,T> ® 

f, 

gi 

.h. 

{g} = I - hjp)/^e (3-72) 

where A* is the master element area. 

F. VALIDATION 

To demonstrate the performance of the finite element formulation developed in the 
present study, the numerical results were compared to the exact solution [Ref 65], and 
existing finite element simulations [Ref 78]. A square smart plate, consisting of a three¬ 
layered (0°/90“/0®) cross-ply laminated plate with the thickness 3 mm was used. A two 
piezoelectric PVDF layers 40 pm each, served as actuator on the top surfece, and 
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a second as sensor on the bottom surfece. The elastic properties that simulate a high 
modulus graphite/epoxy composite are [Ref. 13]. 

Eii= 172.4 GPa (25x10* psi) £ 22 = 6.9 GPa (10* psi); 

Gi 2 =Gi 3 = 3.45 GPa (0.5x10* psi); G23=1.38 GPa (0.2x10*); 

Vi2^13=0.25 

The piezoelectric PVDF layers properties are [Ref 73]; 

Dielectric permittivity 
e 3 i = 0.0460 
e32 = 0.0460 Clm^ 
e33 = 0.0000 C/m'^^ 

Dielectricity 

4 = 0.1062 X lO'** F/m 
4=0.1062x10*’ F/m 
4 = 0.1062 X 10*’ F/m 

Poisson ratio v = 0.29 

Mass density p = 0.1800 x 10^ kg/m ^ 

Modulus of elasticity E=2 x 10 ^ N/m ^ 

The mechanical loading and the electrical potential distribution is described by: 
q = qQ sin(mr /a) sin(;^ / b) 


and 


= V sin(m: / a)sin(; 9 ; / b); 



where is the intensity of load per unit area (N/m \ and F is the anq)litude of O in 
volts. The deflection is normalized as : 


w 


lOOEr 


w 


where; 

Ej- is the transverse Young's modulxis of the graphite/epoxy layers 
A is the span to the thickness ratio 
h is the structure thickness 
=10N/m'^^. 


Figures (3.3), and (3.4 ) show the bending deflections verses length to thickness ratios 
for a smart plate subjected to a double sinusoidal distribution for both mechanical and 
electrical loads, respectively. Figure (3.5) shows a normalized central deflection for a 
sinqjly supported smart plate subjected to a different applied voltage values A three 
dimensional plot of the finite element grid points deflections (81 nodes ) under the same 
loads is shown in Figure (3.6). 

From the results we can conclude that a new finite element model is developed using 
a single hjgher order shear deformation theory to analyze a smart con^she plate. The 
numerical results are conpared with the exact solution [Ref 65], and existing finite 
element simulations [Ref 78]. This result validates the present finite element model. The 
error of the model comes in to the picture at small span to thickness ratios (3 and/or 4) but 
it decreases dramatically as the span to the thickness ratio increases and it attains nearly 
the exact when the ratio exceeds one hundred. The plate deflections can be increased by 
increasing the applied voltage to the actuators. 
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Normalized central deflection 


Central deflection of simply supported smart composite plate 



Length to thickness ratio 

Exact/ [Ref. 65]; Reference a, FE/[Ref 78] 

Figure 3.3: Bending deflection vs. length to thickness ratio for a plate 
subjected to a double sinusoidal mechanical load. 
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Normalized central deflection 


Central deflection of simply supported smart composite plate 



Length to thickness ratio 


Exact/ [Ref. 65]; Reference a, FE/[Ref 78] 


Figure3.4: Bending deflection vs. length to thickness ratio for a plate 

subjected to double sinusoidal electrical and mechanical loads. 
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Normalized central deflection 


Central deflection 



Figure 3.5: Normalized central deflection vs. applied voltage for a plate 
subjected to a uniformly distributed electrical load. 
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Deflection W 


Grid point deflection of simply supported plate(PVDF,(0/90/0),PVDF) 



Figure 3.6: Grid point deflection for simply supported plate subjected to 
a double sinusoidal electrical and mechanical loads. 
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VI. SHAPE CONTROL AND OPTIMIZATION 


A. INTRODUCTION 


Future technologies ai active changes in the structure shape such as airfoils 
spacecraft antenna, turbine blades, etc. justify and necessitate the study of using 
piezoelectric actuators to control the shape of structure deformed quasi-statically. In the 
previous chapter, a finite element model was used to calculate the change in the shape of 
a composite plate t^dien the voltage is applied to the piezoelectric actuators. 

In this part of the analysis, a new model is developed to determine voltage applied to 
the piezoelectric actuator to achieve minimum error between the desired shape and the 
actual shape with prescribed actuator(s) placements. The mathematical model is developed 
based on a kinematical relation between a point in the actual sh^ and a corresponding 
point in the desired sh£q>e for each element of the finite element grid. The sj^em of 
equations for the error function is formulated based on a finite element technique. The 
optimization algorithm is applied to the error (objective) functions through a finite element 
analysis program called OPTSHP. The code performs the analysis using a finite element 
technique and determines the plate deformation due to the ^plied voltage to each 
actuator. It is also able to call a Matlab optimization function to confute the minimiun 
voltages applied to each actuator. 

B. PROBLEM STATMENT 

Consider the problem of a plate with known original shape and actuator configuration 
where the desired shape is specified. Thus it is desired to find actuator voltages to achieve 
the prescribed shape, which minimizes the error between the desired shape and the actual 
shape, for a prescribed actuator(s) position, to achieve a minimum error (objective) 
function. 
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Piezoelectric patch 



Figure 4,1 : Zo coordinate shows a point location on element of the fimte 
element grid 


The analysis is based on small deformation theory. Therefore the specified desired 
shape must be within the region of ‘ small deformation’ firom the original shape. Also, 
the analysis uses the shape change of a reference surfece. Thus within the approximations 
used , the shapes of either the ‘top’ or the ‘bottom’ surfece of the element can be taken 
as the reference surfece. 

The shape of each element of the finite element grid is described via a suitably chosen 
reference surfiice. The shape of this reference surfece is defined by a single -coordinate 
of every point element sur&ce. The Zq- coordinate is perpendicular to the x-y plane as 
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Actual Sur&ce 



shown in Figure (4.1). The shape can be specified a polynomial function in x and y 
such as: 


^0=^0+ ^1^ + + ^3^ + ^4^^ + 

(length units) 

(4.1) 

Also the desired shape can be specified by the same way; 

(length units) 

(4.2) 


^^ere ai and bi are constants. 

C. ERROR FUNCTION A 


Assume that a given set of voltages is applied to a specified number of actuators 
which result in a change in the shape of the plate and give a new configuration of the 
structure, (dashed line in Figure 4.2.) which represent the actual shape. An error function 
will be introduced which includes the difference between the actual shape and the desired 
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one [Ref. 90]. 

The error function will take the form: 

A = l^dxdy (length'’) (4.3) 

where; 

A is proportional to the distance between the reference surfece of the desired shape 
and the actual reference surfiice 
Q, is the element domain. 

A represents the difference between the z-coordinate on the reference surfece of the 
desired shape and the z-coordinate of a corre^nding point on the actual reference 
surfece. 

The error function can also be the sum of the error square at each element 

A = (4.4) 

1=1 

where; A represents the difference between the z-coordinate of a point on the desired 
shape and the z-coordinate of a corresponding point on the actual surfece for each element 
of the finite element grid. 

Thus the objective is to make the error function A as small as possible such that the 
desired surfece and the actual surfece are nearly identical (Le.,yl 0) for a prescribed 
actuator(s) position over the finite element grid, with optimal applied voltages for these 
prescribed optimal placements. 

From Figure 4.3 the distance A is defined as; 

A = (Zo + &o) - i^des + ) (length units) (4.5) 

To evaluate each term in this equation, we consider a nodal point A on the original 
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Actual Surfece 



Figure 4.3: Difference between the actual surfece and the desired surfiice in the x*z plane 
‘TromRef. [91].” 
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reference surfece at (xo,yo). Thus Zo is the z-coordinate of the point A on the original 
reference surfece, is the difference between the z-coordinate of point A on the original 
reference surfece and the z-coordinate of the corresponding point A in the actual surfece , 
which can be ejqjressed as: 


&o = 


--^M + -^V + W 

^dx. 4 ^ 


(4.6) 


where; 

«,v and w are the displacements of the node A. 

w is the displacement tangential to the reference surfece in the axis x-z plane; 

V is the displacement tangential to the reference surfece in the off-axis y-z plane; 
w is the deflection normal to the reference surfece; 

z^ is the z-coordinate of point B on the desired reference surfece at (xo,yo). 


The point on the desired reference surfece corresponding to point A on the actual 
reference surfece is point A located at (Xq + ^o)- ^ difference between 

the z-coordinates of points A and B. 

It be approximated by the first two terms of a Taylor series expansion. 










^0 


(4.7) 


The distances Sx^ and are the changes in the x- and y-coordinates between points A 
andB. 




(4.8) 
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The follow ing par ame ters are introduced which depend on the geometry of the original 
and desired sur&ce, 



By using this parameter the objective function will take the form: 


^ = I , [(^0 - 


(4.10) 


D. FINITE ELEMENT FORMULATION 

A bilinear isoparametric rectangular element is used to describe the shape of the 
element. The element coordinates, x and y are interpolated using a linear shape function 
Nf, which has the form: 

JV-1/4(1 +^4X1+ »7^31) (4.11) 


where; 


Ni is the linear interpolation function. 



^,7/ are the master element coordinates. 
x,y are the actual element coordinates. 

The actual element coordinates x and y can be transformed to the master element by 
the game way as in the finite element analysis. The mathematical relations for the element 
transformation will have the form: 

x = Xc + 4^l2; (4.12) 

y = y, + T}b/2; 

and 

4 = 2faix-x^); 
rj = 2lb(y-y,) 

By substituting equation (4.12) in equations (4.1) and (4.2), the actual and the desired 
shape equations will have the form; 

+ fl,(x, + ^a/2) + a^iyc + >7^2) + «3(^c + + nb/2) (4.13) 

+ a^ix, + ^al2f + a^(y, + vb/2f 
and 

2des=bo+ ^(^c + ^«/2) + b2(yc + nb/2) + bjix^ + ^a/2)(y^ + Tjb/2) (4.14) 
+ b^ (x, + ^al 2 f + b^iy, + 76/2)^ 

(Zq could be zero, for a flat plate, or a first order polynomial) 

By the chain rule of partial diflferentiation, 

&o &„dr] 

3c 3^ & 3ri Sc 

&Q _ 3tj 

dy ^ dy 3rj ^ 
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and 


S^des 

dc dc dt) dc 

&des _^des^ ^de,^ 

^ ^ dn ^ 


(4.16) 


By using equations (4.16) and (4.17), the three geometric parameters qi, ^ 2 , and qs, 
which are e}q)ressed in equation (4.9), can be conq>uted. The objective function equation 
(4.10) takes the form: 


-2<fe(^»7)) + »ll(^.»7)w + »3!!(^.7 )m + (417) 

The electric excitation may be fectored out of displacements u,v, and w in eqiiation 
(4.17) by introducing the vectors j7,vand vv, >^Wch can be defined as follows. The 
degrees of Ifreedom of the system can be determined fi’om the equations of motions (3.57) 
and (3.64) for both cases of actuator and actuator and sensor that are given in chapter 
three in short form as: 

{?}=[r]''{f) (4.18) 

where; 

{?} isthearrayofthedegreesoffteedom andWjj 

[tr]- is the inverse of the structure stifl&iess matrix (given in equations (3.58) and 
equation (3.66)). 

{F*} is the applied electric and mechanical loads vector, in the case of the actuator 
only which is described in chapter three as; 

{r) = {f)+[ji:^][/;;„,]'‘{G) (4.i9) 
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in case of applied electric load only. 

{f"} =[F^iF~.r(G) 

where {G} is the global electric actuating force in which for one element g is defined as; 

{«) = J[ 

and the distributed applied charge density is described as; 




(4.22) 


which is function of the an:q)litude of the electric potential distribution (F in volt). 

The displacement «, v,and w are defined via the interpolation shape function as: 
u = 

/=1 

i=l 

4 

f=l 

«;,Vi,andWj are the displacements at the nodal point in the x, y, and z direction 
respectively, which can be picked out jfrom the deformation array { 9 } expressed in 
equation (4.18) such as: 

«,=J^{F-} (4.24) 

V, = ^{F') 

<*. = f,|f-} 
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where; 

is the flexibility infliience coefficient rows corresponding to the 

displacements ^,v,,andH’^ in the global matrix inverse [a^] 

By substituting u,, v, and w^ equation (4.23) takes the form: 

<4.25) 

<=1 /=1 

i=l 

/=i 

>^iere; 

fj is the shape function used for the deflection w at node j, 

Ni is the shape function used for the displacement u, and v at node j, and w,v,and w 

are the displacements in the three directions, \^ch are confuted fi'om the finite element 
analysis of the plate. 

These displacements are a function of the voltage ^plied to each actuator. To fector 
out the voltage in the objective function, u,v and w, values are used instead of 
u,v, and w in equation (4.17) and one can obtains; 

■^ = |j [(^o(^»'7)-Zdo(^»7)) + + 7ii4,rj)u(y) + 7j(^,7)v(F)f 

(4.26) 

which is a function of the an^litude of the electric potential distribution (F in volt). 
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E. NUMERICAL INTEGRATION 

Equation (4.26) can be integrated numerically using Gauss-Legendre quadrature 
(see Appendix A). A full integration technique of a 3 x 3 Gauss point is used for each 
element. Thus, the error function for the one element can be written as; 

(4.27) 

where j is the determinant of the Jacobian matrix. 

The structure objective function is defined as the summation of the error for each 
element such as: 

i=l 

subjected to 

Lower limit < V < Upper limit 

where m is the number of elements of the finite element grid, and V is the an^Utude of 
the electric potential distribution applied at each actuator (in voltage). 

F. OPTIMIZATION ALGORITHM 

A Matlab function / min and / mins are used to find the minimum of a scalar 
function. The / min function is used to minimize a function of one variable on a fixed 
interval. The problem is mathematically stated as : 


88 




(4.29) 


minimize fix) 

X 

subjected to: 

JC, < X < JCj 

where fix) and x are scalars. 

The /min function is used to optimize the objective function in the a case of constant 
e^lied voltage to each actuator (F is scalar). ' 

The function / mins is used to find the mmimum of a scalar function of several 
variables starting firom an initial estimate, which is generally referred to as unconstrained 
nonlinear optimization, and is mathematically stated as: 

minimize /(x) (4.30) 

X 

x^diere x is a matrix or vector and / is a scalar function. 

The function / mins is used when dififerent values of the voltage are applied to the 
piezoelectric actuators (F is a vector) . / mins uses the simplex search algorithm of 
Nelder and Mead'°^, for minimization of a function of n variables, which depend on the 
comparison of function values at the (n+1) vertices of a general simplex, followed by the 
replacement of the vertex with the highest value by another point. The simplex adapts 
itself to local landscape and contracts to the final minimum. The method is shown to be 
effective and computationally compact. A procedure is given for the estimation of the 
Hessian matrix in the neighbourhood of the minimtun needed in statistical estimation 
problems. 
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G. SOLUTION PROCEDURE 


The flow chart of the solution procedure to confute the optimal voltages applied to 
the actuators and to minimiz e the objective (error) function is shown in Figure 4.4. The 
procedure is as follows: 

1 . Iiqjutthe structural data, such as material properties for both graphite epoxy 
and piezoelectric materials, plate dimensions, actuators dimensions, boundary 
conditions and the t^plied voltages to the actuator. 

2. Assume an initial guess for the applied voltages to each actuator at prescribed 
selected positions. 

3. Call the optimization tdgorithm using aMatlab function /minj (in case of 
different voltages applied to different actuators) or / min (in case of same voltage 
applied to aU actuators). 

4 . Compute the plate deformation due to the current applied voltage, using the finite 
element analysis. 

5. Evaluate the objective fiinction and applied voltages. 

6 . Check the convergence of the objective function and voltages. If they converge, 
the program will terminate. 

7. Otherwise, update the current applied voltages to the actuator(s) and then go to 
step three and repeat the procedure. 

H. VALIDATION 

As illustrated in the flow chart. Figure 4.4, the optimal location of the actuator(s) 
position to get a minimum error function at a specified applied voltage is determined. 
The OPTSHP code for the optimization algorithm has been tested. 
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Figure 4.4: Flow chart of solution procedure 
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Figure 4.5 Element numbers for different actuator positions for a sin:5)ly supported 
square con^site plate (a= b). 


A square fiber-reinforced con^sHe plate with three layers (0**/90*’/0®) and one 
piezoelectric layer patch on the top surfiice and another one on the bottom surfece was 
used. The length to thickness ratio equal fifty, and the materials properties for 
graphite/epoxy are as follows [Ref 13]. 

Eii= 172.4 GPa (25x10^ psi) E22=6.9 GPa (10^ psi); 

Gi2=Gi3 =3.45 GPa (0.5x10® psi); G23=1.38 GPa (0.2x10®); 

Vi2=Vi3=0.25 


The piezoelectric PVDF layers properties are [Ref 73]: 
Dielectric permittivity 

esi = 0.0460 C/m''^ 
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632 = 0.0460 C/m^ 

633 = 0.0000 C / m ^ 

Dfelecstricity 

=0.1062 X 10-® F/m 
4=0.1062x 10-’ F/m 
4=0.1062x10'’ F/m 

Poisson ratio v = 0.29 

Mass donsity p = 0.1800 x 10^ kg/m^^ 

Modulus of 6lasticity E=2 x lO"^ N/m^^ 

Th6 mochanical loading is takon to bo zero, and tho oloctrical potontkl distribution is 
assumed to be uniform^ distributed with anqilitude V. The deflection is normalized as : 

_ 100£r 

H- =- —W 

where; 

Ef is the transverse Young's modulus of the graphite/epoxy layers 
X is the span to the thickness ratio (equal fifty) 
h is the structure thickness 

A nine element model was used for different positions of the actuators as shown in 
Figure 4.5. The original sh^ is chosen as a flat plate {z^ = 0) and the desired shape is 

selected as: = 1 x 10“’- 42 x 10“’®+ 6 x 10“V^ • 


93 


Exan^le 1. 

In this case the error between the actual sh^ and the desired shape due to a specified 
amount of the applied voltage to the actuator is computed by placing the actuator(s) at 
different positions on the finite element grid (Figure 4.5 ). The error function values are 
given in Table 4.1 for an ^plied voltage equal to 100 volts for each actuator position. 

Exan^le 2. 

The predicted voltages and the grid points transverse deflections, for both the actual 
and the desired transverse deflections, (in-plane deformations are not included) are 
con 5 )uted for two cases. Figure (4.6) shows the first case, when the actuator is placed at 
element number five. The second case, vriien two actuators are used at elements two and 
eight, is shown in Figure 4.7. 

Exan:5)le 3. 

The optimal value of voltage applied at the actuator(s) to obtain minimum error 
between the actual shap e and the desired shape (minimum error fiinction) can be 
enmpiited as shown in the fiow chart. Figure 4.4, using the Matlab f min function with 

certain upper and lower limits of the applied voltage, such as -100< V < 100. The optimal 
voltage jqiplied to the actuator for different positions on the grid and the corresponding 
minimum value of the error function are given in Table 4.2. 

The reported values were the global minimums. This was acconqjlished by applying 
different initial values of voltages (Le. -500 to 500 Volt.) to the actuator and determining 
the n ume rical value of the objective function at these different voltages. The m i n i m u m 
value of the error function found was identical to the value computed by the optimization 
algorithm. 
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Table 4.1: Error function values at different actuator positions for applied voltage equal 


100 vok. 


Actuator Position 

Error Function Value 

Element # 1 

f = 5.87619804 e-18 

Element #2 

f = 5.05453907 e-18 

Element #3 

f = 4.44646994 e-18 

Element #4 

f = IMOimS e-18 

Element #5 

f = 9.30776749 e-18 

Element #6 

f= 4.28947367 e-18 

Element #7 

f = 5.86952426 e-18 

Element #8 

f = 5.04220635 e-18 

Element #9 

f = 4.44063580 e-18 


• The applied ventage at each actuator is 100 Volts. 


Example 4. 

In the case of four actuators used at a time, the voltages and error function values are 
given in Table 4.3. Table 4.4 shows the results for two actuators used at a time. The 
results in both cases were obtained by using the Matlab function / min .v with 
unconstrained value on the ^plied voltage. 

An optimization algorithm is applied to get minim um applied voltages to the actuators to 
minimize the error function. Tables 4.3 and 4.4 show that location of the actuators play an 
in^rtant role in minimizing the error. 
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Figure 4.7: Actual and desired transverse deflection at the grid point 
for actuators at elements two and eight. 
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Table 4.2: Optimal iq)plied voltage and error function at different actuator positions. 


Actuator Position 

Minimum Applied 

Voltage (Voh.) 

Minimum Error 

Function 

Element # 1 

V = 78.226 

f=5.83495108e-18 

Element #2 

V = 70.944 

f=4.78995246e-18 

Element # 3 

V = 99.999 

f=4.44647064e-18 

Element #4 

V = 34.161 

f=5.88349351e-18 

Element # 5 

V = 37.346 

f=4.74690856e-18 

Element #6 

V = 75.062 

f=4.03165923e-18 

Element #7 

V = 78.611 

f=5.82972138e-18 

Element # 8 

V = 71.141 

f= 4.78117803e-18 

Element # 9 

V = 99.999 

f=4.44063651e-18 

Table 4.3: Optimal applied voltages and error fimctions for the case of four actuators used 

at a time. 

Actuators Positions 

Minimum Applied 

Voltages (Volt.) 

Minimum Error 

Function 

Elements #' 1,3, 7 & 9 

Vi=-21.060 

Vs = 98.104 

V 7 =-16.473 

Vs = 101.482 

f= 4.0531138e-18 

Elements 2,4, 8 & 6 

V 2 = 18.404 

V 4 = -29.429 

V 6 = 65.262 

V 8 = 22.349 

f= 4.01531283e-l 8 
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Table 4.4: Optimal applied voltages and error functions for the case of pair of actuator 
used at a time. 


Actuators Positions 

Minimum Applied 

Voltages (Voh.) 

Minimum Error 

Function 

Elements #* 1 & 7 

Vi =41.612 

V7 = 47.265 

f= 5.78007896e-18 

Elements #’ 2 & 8 

Vz* 35.943 

Vg = 39.934 

f=4.74678213e-18 

Elements #* 3 & 9 

V3 = 85.710 

Vs = 90.547 

f=4.0583368e-18 

Elements #* 1 & 3 

V,=-26.368 

Vs = 175.055 

f= 4.10114654e-18 

Elements #* 4 & 6 

V4 =-8.231 

Vs = 79.660 

f=4.01544062e-18 

Elements #* 7 & 9 

V, = 26.0895 

Vs = 175.193 

f=4.09288468e-18 

Elements #* 1 & 9 

V, =-17.007 

Vs =170.053 

f=4.10970645e-18 

Elements #* 3 & 7 

Vs = 169.712 

V7 =-16.138 

f=4.12087738e-18 
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V. CONCLUSIONS 


A finite elen^nt model was developed to analyze the behavior of a smart con^site 
plate. This intelligent plate was con^sed of distribxited sensor and actuator layers made 
of PVDF or PZT and a laminate of graphite/epoxy layers. In the analysis, the piezoelectric 
layer was treated as a normal lamina. A siipple higher order shear deformation theory 
with Hamilton's princh)le was used to formulate the equations of motion. A four node, 
bilinear, isoparametric, rectangular element, with seven degrees of fieedom at each node 
was developed. The electric potential was treated as a generalized electric coordinate like 
the generalized displacement coordinates at the mid-plane of the actuator and sensor 
layers. A Matlab code ‘CMPZ’ was developed to analyze the structure. 

Several conclusions can be made fi'om the results. The present finite element method 
can be used for static and dynamic analyses for both thick and thin con^she structures 
with distributed piezoelectric sensors and actuators. The numerical results generated by 
the developed code agree very well with the exact solution and other published finite 
element solutions. This validates the finite element model The method developed is much 
sinqjler to formulate and more efficient than models based on hexahedral or tetrahedral 
solid elements and/or three dimensional brick elements. The error of the method 
increased for small span to thickness ratios and decreased dramatically for higher span to 
thickness ratios. A Hermite cubic interpolation function was used to approximate the 
transverse deflections, however, the method does not suffer fi-om the shear correction, 
\^bich is problematic in the first order shear deformation theory. The developed 
displacement model includes the parabolic distribution of the transverse shear stresses and 
the non-linearity of in-plane di^lacements across the thickness. This is an advantage over 
the classical laminated plate theory, which neglects the effects of transverse shear stresses. 
The number of degrees of fi-eedom of the element used in the present method is one third 
the number of degrees of fi'eedom of the element used in the model developed by the 
higher order shear deformation theory, which of course saves con^utational time. 
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For shape control and optimization, an optimization algorithm based on a finite 
element method was developed. Kinematical relations were used to formulate the 
objective (error) function as the summation of the mean square error between a point in 
the actual sur&ce and a corresponding point in the desired surfece integrated for each 
element of the finite element gnd. This approach is an in^rovement over the method 
where surfiice error is determined only at the nodal points. 

The optimal voltage applied at each actuator to achieve a specified shape with 
rriinimiim eiTor between the actual shape and the desired shape was conqjuted using an 
optimization code. The Matlab code ‘OPTSHP’ is developed in conjunction with Matlab 
functions which gave very satisfectory results. The code is able to determine the value of 
the error function for different actuator(s) positions for a specified amount of the applied 
voltage to the actuator(s). The model is applied for both cases, either constant voltage 
applied to the actuator(s) or different voltages applied to several actuators. The 
procedure proposed in this work was inqilemented for a composite plate with any 
demonstrated number of piezoelectric actuators. It was demonstrated firom the examples 
that the desired shape was made to chosen match the actual sh^ie, 'wdiich satisfies small 
deformation theory. The procedure was tested for more general layouts of actuators and 
the optimization algorithm was applied to get optimal applied voltages to the actuators to 
minimize the error function. 

Future work is proposed in the following areas; 

1. Generalization of the finite element model to non-rectangular elements and shell 
elements. 

2. Optimization of the actuator placements. 

3. Closed loop shape control to compensate for varying external mechanical and thermal 
loads. 
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APPENDEXA. NUMERICAL INTEGRATION 


For the isoparametric element considered in this work, an exact evaluation of the 
integration for mass, stffiiess and the equivalent nodal loads was not always possible 
because of the algebraic conqjlexity of the differential equations representing the element 
used. Since the interpolation functions are easily derivable for a rectangular element, and 
it is easier to evaluate integrals over rectangular geometry, we transform the finite element 
integral statements, defined over quadrilaterals, to a rectangle. In such a case, it is 
necessary to use some numerical integration technique. One of the most accurate and 
convenient metlK>ds is Gaussian quadrature, which involves approximation of the 
integrand by a po^omial of sufficient degree, because the integral of a polynomial can be 
evaluated exactly [Ref. 93]. If an integral; 

P=J['F(x)d&r (Al) 

with function F(x) is approximated by a polynomial; 

N 

(A-2) 

\\1iere FJ denotes the value of F(x) at the I* point of the interval (xi,X 2 ) and are 
polynomials of degreeA-1. The representation can be viewed as the finite element 
interpolation of F(x) , where F] is the value of the function at the I* node. For our 
rectangular element the integral will take the form; 

M N 

[ F(x,y)dx<iy=l ^ 

" ' /=1 y=l 

(A.3) 
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where; 

is the actual element domain with its x, and y coordinates ( Fig. 2.3) 
is a linear master element with its r and s coordinates (Fig.2.3) 

(^, rf) denote the Gauss points in the ^ and 17 directions. (Fig 2.3) 

R, and Rj denote the corresponding Gauss weights (see Table. A.1) 

A/and N denote the number of quadrature points. (in most, cases, the interpolation 
functions are of the same degree in both % and ti direction, Le. M — N. 

I J1 is the determinant of the Jacobian matrix J, in ^^ch ,| > 0 every\^ere in the 

element ^2^. 

Geometrically, the Jacobian represents the ratio of an area element in the actual 
element to the corresponding area element in the master element: 

dA = dxcfy = \j\d^7j (A.4) 

The number of Gauss points is based on the formula; 

N = mteger(lf2(p + 1)) (A.5) 

which means that the smallest integer greater than 1/2 (p + 1 ); >^iiere p is the degree of the 
polynomial which may be integrated exactly enqjloying equation (A.5). 

Table A.1 gives the integration order, weighting fector and the location of the Gauss 
points for linear, quadratic, and cubic elements [Ref 94]. The maximum degree of the 
polynomial refers to the degree of the highest polynomial in 4 and t) that is present in the 
integral F{^,Tf) of equation (A.3). 
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Table A.1: 

Weights and Gauss points for the Gauss- Legendre Quadrature. 


N 


1 

1 

0.0 

2.0 

2 

0.5773502692 

1.0 

3 

0.7745966692 

0.5555555556 


0.0 

0.8888888889 

4 

0.8611363116 

0.3478548451 


0.3399810436 

0.6521451549 

5 

0.9061798459 

0.2369268851 


0.5384693101 

0.4786286705 


0.0 

0.5688888889 

6 

0.9324695142 

0.1713244924 


0.6612093865 

0.3607615730 


0.2386191861 

0.4679139346 

7 

0.9491079123 

0.1294849662 


0.7415311856 

0.2797053915 


0.4058451514 

0.3818300505 


0.0 

0.4179591837 

8 

0.9602898565 

0.1012285363 


0.7966664774 

0.2223810345 


0.5255324099 

0.3137066459 


0.1834346425 

0.3626837834 


N=2, Linear, N=3, Quadratic; N»4 Cubic, etc. 

N X N *the order of integration 
% - Location of mtegraticm points in master element. 
R» Weighting fector. 
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APPENDEXB. COMPUTER PROGRAMS 


Two Matlab codes were developed, 'COMPZ' and 'OPTSHP'. The 'COMPZ' code 
is a finit e element analysis program for a smart conqtosite plate. The program is able to 
co mp ute the Static and the dynamic response of the structure, bending deformation, free 
vibration, and stress analysis subjected to both mechanical and electric load for various 
boundary conditions. The program is able to solve a composite plate with and without 
piezoelectric layers. The composite laminates could be graphite epoxy with reinforced 
fiber or ahimimim layers. The piezoelectric elements can be segmented or continuous 
elements and can be either surfoce bonded or embedded in the laminated plate. The 
procedure to run the program is as follows: 

1. Mesh generation: 

Input the number of point in x direction Nx 

Input the number of point in y direction Ny 

Input the plate length in X direction Lx 

Ii^ut the plate length in y direction Ly 

Input the plate thickness in z direction Lthk 

Input the number of patches their lengths in x direction (if any) 

Input the number of patches their lengths in y direction (if any) 

The program will generate the mesh automatically and compute the elements numbers, 
nodes coordinates, and the elements conductivity. 

2. Structures Data: 

Input type of the structures, singly supported or cantilever plate. 

Input materials constants for the structures core (graphite epoxy or aluminum) and for 
the piezoelectric material. E), E:, etc. 
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3. Applied Loads: 

Input the intensity of the load per unite area q©. 

Ii^)ut the amplitude of the {q)plied voltage V. 

4. Method of analysis: 

Input quarter plate analysis. 

Input full plate analysis. 

Con:q>osite plate analysis. 

Smart con^site plate anafysis, with piezoelectric materials. 
Complete cover the sur&ce, or 
Uniform patches cover the sur&ce 

5. Resuh 

Plate deflection 

Natural frequency and mode shape 
Stress value for each layer 


The "OPTSHP" code is an optimization program using the finite element technique in 
which the optimization algorithm is applied to the objective function by using a Matlab 
function / mins through finite element analysis subroutines. The code is able to compute 
the plate deformation due to the applied voltage and/or mechanical load and determine the 
optimal placement of actuator. Also the code is able to conpute the minimum amoimt of 
the applied voltages to each actuator by using a Matlab function /mins and /min. 

The steps to run the code are: 

1. Mesh generation: 
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Input the number of point in x direction Nx 
Ii^)ut the number of point in y direction Ny 
li^ut the plate length in X direction Lx 
Input the plate length in y direction Ly 
li^ut the plate thickness in z direction Lthk 
Input the number of patches 
Input the patches element numbers 

2. Structures Data: 

It 5 )ut type of the structures, singly supported or cantilever plate. 

Input materials constants for the structures core (graphite epoxy or aluminum) and for 
the piezoelectric material. Ei, E 2 , etc. 

3. Applied Loads: 

Input the intensity of the load per unite area qo (if any) 

Ii^ut the initial value of the applied voltage V . 

4. Method of analysis: 

Input full plate analysis. 

Smart composite plate analysis, with piezoelectric materials. 

Uniform patches cover the surfece 

5. Determine the patch location: 

Input initial guess for patch position and conq)ute the objective function 

The program will do finite element analysis and compute the plate deflection due 

to the applied voltage at this place. 

Input the second guess at other location of the finite element grid. 
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Repeat this dieck and corapaie the values of the objective function for each case and 
pick the minimum value which represent the optimal place to get minimum error 
between the actual and the desired shape. 

6. Determine the optimal amount voltages {q>plied to each actuator: 

Input initial value of the applied voltages for the each actuator at the selected position. 
Call the optimization algorithm using the Matlab function f mins and 
/ min, they will perform a finite element analysis for this value of the voltages. 

The code will update the voltage values and repeat the analysis again. 

The code will check the global minimum of the objective function and will stopped at 
the global minimum with minimum amount of voltages s^plied to each actuators. 
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